#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>
#include <fstream>
#include "QuadraticBSplines.h"
#include "CubicBSplines.h"

double f(double x) {
    return 1.0 / (1.0 + pow(x, 2));
}

int main() {
    double m_1, m_n;
    double epsilon = 1e-6; // 更小的epsilon值以提高精度
    m_1 = (f(-6 + epsilon) - f(-6)) / epsilon;
    m_n = (f(5) - f(5 - epsilon)) / epsilon;
    std::vector<double> x1, f_value1;
    double f_1 = f(-11/2.0);
    for (int i = 0; i < 11; i++) {
        x1.push_back(-6.0 + i);
        f_value1.push_back(f(x1.back()));
    }

    std::vector<double> x2, f_value2;
    double f_n = f(11 /2.0);
    for (int i = 0; i < 11; i++) { // 从0开始以匹配x1的索引
        x2.push_back(-9 / 2.0 + i); // 修正x2的计算
        f_value2.push_back(f(x2.back()));
    }

    std::string dataFilename = "spline_data_cubic.txt";
    CubicBSpline_complete spline_complete(x1, f_value1, m_1, m_n);
    std::ofstream outFile(dataFilename);
    int numDensePoints = 500;
    for (int k = 0; k < numDensePoints; k++) {
        double xi = -6.0 + 11.0 * k / (numDensePoints - 1);
        double y_true = f(xi);
        double y_complete = spline_complete.B_value(xi);
        outFile << xi << " " << y_true << " " << y_complete << "\n";
    }
    outFile.close();

    std::string scriptFilename = "plot_splines_cubic.gp";
    std::ofstream scriptFile(scriptFilename);
    scriptFile << "set terminal pngcairo enhanced color size 1200,800\n";
    scriptFile << "set output 'spline_comparison_cubic.png'\n";
    scriptFile << "set title 'Comparison of Splines for Thm3.57 '\n";
    scriptFile << "set xlabel 'x'\n";
    scriptFile << "set ylabel 'y'\n";
    scriptFile << "set key outside\n";
    scriptFile << "plot '" << dataFilename << "' using 1:2 with lines title 'True Function', \\\n";
    scriptFile << "     '" << dataFilename << "' using 1:3 with lines title 'Complete B Spline'\n";
    scriptFile.close();
    system(("gnuplot " + scriptFilename).c_str());

    std::string dataFilename2 = "spline_data_quadratic.txt";
    QuadraticBSplines spline_quadratic(x2, f_value2, f_1, f_n);
    std::ofstream outFile2(dataFilename2);
    for (int k = 0; k < numDensePoints; k++) {
        double xi = -9 / 2.0 + 9.0 * k / (numDensePoints - 1); // 修正xi的计算以反映二次B样条的定义域
        double y_true = f(xi);
        double y_quadratic = spline_quadratic.B_value(xi); // 使用正确的B样条对象
        outFile2 << xi << " " << y_true << " " << y_quadratic << "\n";
        }
outFile2.close();
    outFile2.close();

    std::string scriptFilename2 = "plot_splines_quadratic.gp";
    std::ofstream scriptFile2(scriptFilename2);
    scriptFile2 << "set terminal pngcairo enhanced color size 1200,800\n";
    scriptFile2 << "set output 'spline_comparison_quadratic.png'\n";
    scriptFile2 << "set title 'Comparison of Splines for Thm3.58 '\n";
    scriptFile2 << "set xlabel 'x'\n";
    scriptFile2 << "set ylabel 'y'\n";
    scriptFile2 << "set key outside\n";
    scriptFile2 << "plot '" << dataFilename2 << "' using 1:2 with lines title 'True Function', \\\n";
    scriptFile2 << "     '" << dataFilename2 << "' using 1:3 with lines title 'quadratic B Spline'\n";
    scriptFile2.close();
    system(("gnuplot " + scriptFilename2).c_str());
    std::cout<<"figure saved as spline_comparison_cubic.png and spline_comparison_quadratic.png"<<endl;
    std::vector<double> points = {-3.5, -3, -0.5, 0, 0.5, 3, 3.5};
    std::vector<double> errors_complete, errors_quadratic;
    for (int j=0;j<points.size();j++) {
        double x = points[j];
        double y_true = f(x);
        double y_complete = spline_complete.B_value(x);
        double y_quadratic = spline_quadratic.B_value(x);
        errors_complete.push_back(fabs(y_complete - y_true));
        errors_quadratic.push_back(fabs(y_quadratic - y_true));
    }

    // 输出插值误差
    std::cout << std::fixed << std::setprecision(6);
    for (size_t i = 0; i < points.size(); ++i) {
        std::cout << "Point: " << points[i] 
                  << ", Cubic Error: " << errors_complete[i] 
                  << ", Quadratic Error: " << errors_quadratic[i] << std::endl;
    }

    return 0;
}